# ============================================================
# Thomas König & Stefan Eschenwecker
# The European Court of Justice and legal European integration
# 
# This script produces the regression model object 
# for the No MS Model (see Table A3 in online appendix).
# ============================================================

# load packages
library(rstan)
library(brms)
library(tidyverse)

# set options for the script
options(mc.cores = parallel::detectCores())
options(scipen = 10)


# load data
CourtData <- read.csv("Court/Data/PreparedCourtData.csv")

# prepare data
# (DV + IVs as ordered factors)
# (remove cases where all actors claimed not applicable)
CourtData <- CourtData %>%
  mutate(
    pos_cjeu_ln = factor(pos_cjeu_ln,
      ordered = T,
      levels = c("PS", "Ambi", "ME")
    ),
    supra_signal = factor(supra_signal,
      ordered = F,
      levels = c( "Uninfo", "Cont", "Weak PS", "Strong PS", "Weak ME", "Strong ME")
    ),
    supra_signal_int = factor(supra_signal_int,
      ordered = T,
      levels = c("No Intensity", "Weak Intensity", "Strong Intensity"))
  ) %>%
  filter(all_not_applic == 0)


# estimate regression model
# (automatically save the object)
OrdProbModel <- brm(
  bf(
    pos_cjeu_ln ~ supra_signal +
      share_judges_std + (1 | iuropa_case_id)
  ) +
    lf(
      disc ~ 0 + mo(supra_signal_int) +
      share_judges_std + (1 | iuropa_case_id),
       cmc = F),
  family = cumulative(link = "probit"),
  prior = c(
    # default alpha = 1 Dirichlet prior for monotonic effect
    prior(normal(0, 1), class = b), # standardized location predictors
    prior(exponential(1 / 0.463), class = sd, group = iuropa_case_id), # 95% interval of 0.25 to 4
    prior(normal(-0.43, 2), class = Intercept, coef = 1), # location is qnorm of 1/3
    prior(normal(0.43, 2), class = Intercept, coef = 2), # location is qnorm of 2/3
    prior(normal(0, log(4) / 2), class = b, dpar = disc), # expect difference by factor of max 4
    prior(exponential(1 / 0.463), class = sd, dpar = disc) # 95% interval of 0.25 to 4
  ),
  init = 0, # avoid rejected starting values
  data = CourtData,
  chains = 4, cores = 4, warmup = 1000, iter = 2500,
  seed = 1408, control = list(adapt_delta = 0.99),
  file = "Court/Models/UnequalPrecision_NoMS",
  backend = "cmdstanr"
)

# print summary
summary(OrdProbModel)

# check for non-convergence
shinystan::launch_shinystan(OrdProbModel)
